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We solve the model of quantum Brownian oscillators linearly coupled to an environment of 
quantum oscillators at finite temperature, with no extra assumptions about the structure of the 
system-environment coupling. Using a compact phase-space formalism, we give a rather quick and 
direct derivation of the master equation and its solutions for general spectral functions and arbitrary 
temperatures. Since our framework is intrinsically nonperturbative, we are able to analyze the 
entanglement dynamics of two oscillators coupled to a common scalar field in previously unexplored 
regimes, such as off resonance and strong coupling. 
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I. INTRODUCTION 
A. NQBM 

In Ref. [1] exact solutions were obtained for the prob- 
lem of one Brownian oscillator linearly coupled to an 
environment of quantum oscillators, or quantum Brow- 
nian motion (QBM). It was briefly mentioned how the 



phase-space formalism therein rather trivially general- 
ized to the problem of N Brownian oscillators with arbi- 
trary spectral-density function, or NQBM. In this work 
we make explicit this generalization and clarify details 
which become more important in multipartite systems. 

An important discovery made in Ref. [ ] was that pre- 
vious derivations of the QBM master equation were in- 
valid for nonlocal dissipation, due to a subtle issue con- 
cerning boundary conditions of integro-differential equa- 
tions. This issue remains when considering multiple sys- 
tem oscillators and, thus, applies also to NQBM treat- 
ments prior to Ref. [ ]. Another important result there 
was the calculation of exact, analytic results for sub- 
Ohmic, Ohmic, and supra-Ohmic couplings. The suit- 
ability of our QBM formalism for such exact calculations 
remains in the NQBM case. However, as a nontrivial 
application of our formalism we study here the case of 
several detectors coupled to a relativistic-field environ- 
ment, where the spectral density is far more complicated 
than a simple power law. 

NQBM was previously considered in Ref. [ . ] , but only 
for couplings such that the environment is effectively cou- 
pled only to one system degree of freedom. Ref. [ ] went 
into great detail for resonant system oscillators, but oth- 
erwise treated the problem generally enough to model 
local oscillators interacting via a relativistic field. The 
most general regime was briefly discussed in Ref. [l] and 
has since been covered by Ref. [4] in some detail. In this 
work we develop the formalism and expressions which al- 
lowed us to obtain the exact analytic solutions in Ref. [ ] . 



B. Entanglement Dynamics 

In the second part of the paper we apply our formal- 
ism to investigate the entanglement dynamics of local 
oscillators in a common scalar- field environment. The 
scalar-field environment provides a source of noise which 
is correlated, not only in time, but also in space. This 
detail has a large effect upon the entanglement dynamics 
of the oscillators. The regime we consider is very general 
and encapsulates those considered in other treatments. 
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The study in Ref. ["] essentially corresponds to the case 
of two local and resonant oscillators at precisely the same 
location within their environment, such that their en- 
vironmental influences were identical. Ref. ['■] carried 
out a perturbative analysis of the entanglement dynam- 
ics of local and resonant oscillators at separate locations 
in a scalar- field environment, so that their environmental 
influences were intricately correlated according to their 
separation. On the other hand, Ref. [4] considered non- 
resonant oscillators but with a generalized Ohmic cou- 
pling rather than a field environment exhibiting nontriv- 
ial spatial correlations. 



C. Organization 



to inserting the entire system-environment interaction in 
the square of the potential: 

L= ^ (X'^MX-X'^CX) (II.3) 
+ i (k^m±- [x + c-igXj'^c [x-Hc-igX]) , 

since this keeps the phenomenological system-system 
couplings from changing. Ref. [A] chose to only renor- 
malize the environment ally- induced self interactions or 
diagonal terms, which are the most cutoff sensitive terms. 
In some ways that choice is more analogous to quantum 
electrodynamics . 



In Sec. II we provide a fairly self-contained derivation 
and analysis of the A^-oscillator solutions. This includes 
solutions to the Langevin equations and master equa- 
tion. Due to the matrix formulation in Ref. [1], most 
of the useful relations therein can be generalized quite 
naturally. In Sec. Ill we develop relations for the system 
evolution and late-time covariance which are especially 
useful in multipartite systems. Finally in Sec. IV we ap- 
ply this formalism to the problem of a pair of local os- 
cillators residing in the same scalar field, but with some 
finite separation between them. 

Throughout the paper we use natural units with h — 
c= 1. 



II. N-QBM 

A. The Lagrangian 

Our model is that of a continuous and linear system 
with finite and countable degrees of freedom, with La- 
grangian Lsys(X, X), bilinearly coupled, via a Lagrangian 
Lint(X,x), to a linear environment with an infinite (and 
possibly continuous) number of degrees of freedom, with 
Lagrangian Lcnv(x,x). 

L = Lsys(X, X) + Lc„v(x, x) + Lint(X, x) + Li.c„(X) , 

(ILl) 

L = i (X'^MX-X'^CX) + ^ (i^mi-x^cx) 

-xTgX + Lre„(X). (IL2) 

We assume that the spring constant matrices C, c as well 
as the mass matrices M, m are real and positive definite, 
and can be considered in general to be symmetric. If nec- 
essary, one can relax the positivity condition and even 
consider time-dependent mass matrices, spring constant 
matrices and system environment coupling matrix g [' ]. 
To ensure that the free and interacting system are simi- 
lar in behavior, we will also include the renormalization 
Lion(X). Our choice of renormalization will be equivalent 



B. The Langevin Equation 

For the linear system there are several formalisms 
which produce the same Langevin Equation. The most 
direct is via integrating out environment degrees of free- 
dom in the Heisenberg equations of motion [7] and then 
considering the symmetrized moments. Another is to 
consider the characteristic curves of the system -I- envi- 
ronment's Fokker-Plank equation [()]. Finally, one can in- 
tegrate out both the environment degrees of freedom and 
the relative system coordinate A = X — X', while leaving 
only the average system coordinate S = (X — X')/2, in 
the double path integral of the reduced system propaga- 
tor in the influence functional formalism [ ] . In general 
(for nonlinear systems) there is no necessary correspon- 
dence between these formalisms and only the first may 
be well defined, but here the Langevin equation is simply 

M±{t) + {C + 5C)X{t) + 2 f dr ti{t, t) X(r) = ^ (<) , 

Jo 

(11.4) 

where 6C is a renormalization of the system interac- 
tion, /i(t, t) is the dissipation kernel, and ^{t) a Gaussian 
stochastic force with vanishing mean and correlation 

(meir))^^i^it,r), (II.5) 

characterized by the real noise kernel i'{t,T). In the 
phase-space representation, ^(t) can be formally inter- 
preted as a Gaussian stochastic process and X(t) as an 
unraveling of the characteristics which evolve the quan- 
tum state. 

In general, for any Gaussian environment (and not only 
the one we model here) the two kernels are constrained 
by the fact that the environment correlator 

a(t,r) = /y(t,T) +i/x(i,T), (II.6) 

must be Hermitian and positive definite ['']. Specifically, 
for a constant bilinear coupling to a thermal reservoir of 
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harmonic oscillators we have stationary kernels 
, , _ 1 sinfwft — rl) _i 

Z U) 



(II.7) 



T 1 coth(:|^) cos(a;[i-T]) i 
iy(t,r) = +g^m^5 ^^^^^^ ^-^ ^m"5g, (Il.i 



2a; 



o;^ = m 2 cm 2 ^ 



(II.9) 

with temperature-independent dissipation due to linear- 
ity in the field coupling. 

For the positive-temperature environment, dissipation 
may also be represented via the positive-definite damping 
kernel defined by the relation fJ-{t,T) = — ^7(t,T), such 
that 

T 1 cos(u;[i — rl) 1 . ^ 

7^,T =+g^m~2 y ^V -2g, 11.10 

and with which the Langevin equation can be expressed 
as 

MX(t) + 2 [ dr j{t, t) X(r) + CX(t) (11.11) 
Jo 

+2 Jit) Xo^m, 

The noise and damping kernels satisfy then the 
fluctuation-dissipation relation (here in the Fourier do- 
main) 

i>{uj) ^ k{Lj)-y{uj) , (11.12) 
k{oj) = hw coth(^^^Y (11.13) 
with the Fourier transform defined 



/» 



+ 00 

dte'"^' 

— 00 



(11.14) 



and where k is the (quantum) FDR kernel. Therefore, the 
problem is completely specified in terms of the damping 
kernel (or equivalently, the spectral-density function). 

Given a stationary damping kernel, the Langevin equa- 
tion can then be expressed in the Laplace domain as 



[s^M + 2sj{s) + C] X(s) = [sMXo + Pq] 



(11.15) 



where P = MX and (Xo,Po) correspond to the initial 
values at i = 0. Formally, the solutions can be easily 
found by inversion: 

X(s) = G(s) [sMXo + Po] + G{s) ^(s) , (IL16) 

G(s) = [s^M + 2sj{s) + C]^\ (11.17) 

Note that since our damping kernel is symmetric, i.e. 
7(t, t) — 'y'^(t, r), the same will be true for the propaga- 
tor G(t,T) and its Laplace transform. It is also useful to 
consider the following representation: 

.1-1 



G(s) = M-2 + 2s A(s) + 
X(s) = M-^(s)M-^ , 



M 



M-2CM~2 



(IL18) 

(IL19) 
(11.20) 



where the eigenvalues of coincide with the squared 
frequencies of the normal modes of the free system. Back 
in the time domain we have 

X(t) = G(t) M Xo + G(t) Po + (G * m) , (11.21) 

with * denoting the Laplace convolution, defined as 

(A*B)(t)^ [ dTA(t-T)B(T). (11.22) 
Jo 

At this point it is worth pointing out that prior to [1, 
4], non-perturbative iV-oscillator master equations have 
been limited to the parameter regimes in which the A'^- 
oscillator problem was equivalent to independent single- 
oscillator problems. More specifically, either the coupling 
was chosen such that 7 c>c 1 [ ] or the the resonance 
regime with C oc 1 was considered [-], or both [5]. It 
is more or less obvious from the Langevin equation that 
these examples, where 7 and C commute, are all trivial. 

C. Phase Space: The Solutions 

The compact and powerful formalism developed in 
Ref. [1] and which allowed a quick derivation of the mas- 
ter equation and its general solution, also makes straight- 
forward the generalization to N system oscillators. First, 
one introduces the following phase-space representation 



N(i,T) 

T(t,r) 



Z^= (X,P) , 
3^^(0,0 , 




u(t,T) 


j(t,T) 



(11.23) 
(11.24) 

(11.25) 
(11.26) 



The solutions of the Langevin equation can then be ex- 
pressed as 



Z(t) = ^(t)Zo + (^*S)(t) 

G(t)M G(t) 
MG(t)M MG(t) 



(11.27) 
(11.28) 



As shown in Ref. [ ] , the reduced Wigner function of the 
system can be conveniently written as a double average 
over the stochastic source and the initial conditions: 



W(Z,t) = 



(11.29) 



where Z(t) is a solution of the Langevin equation as given 
by (11.27). It is then convenient to consider the Fourier 
transform from phase-space coordinates Z to the conju- 
gate Fourier-domain coordinates K 



W(i,K) 



d^Xd'^Pe~"'^''^W(t,Z), (11.30) 
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which corresponds to the characteristic function from 
which the Weyl-ordered correlation functions can be di- 
rectly obtained by functionally differentiating with re- 
spect to K. Following the same procedure as in Sec. 2.3 
of Ref. [ ] , one can immediately get 

W{t, K) = W(0, *'^(<)K) g-iK^-^^WK ^ (11.31) 

where crxit) is the thermal covariance matrix 

arit)^ [ dr [ dr' *(<-T)N(r, r')*'^(i-r') . (11.32) 
Jo Jo 

Since the Wigner function is a pseudo-distribution (it 
is real, normalized, but not necessarily everywhere- 
positive) its Fourier transform is a characteristic (or mo- 
ment generating) function. A straightforward interpre- 
tation of the solution immediately follows: the initial 
cumulants of the system experience damped oscillations 
due to the homogeneous propagator $(t) while a thermal 
covariance arises, smearing out details of the dissipating 
initial state. 

The thermal covariance can also be broken down into 
the block-matrix correlations 

(T^^{t) ^ I dr [ dr' G(t~T)u(T,T')G^{t-T') , 
Jo Jo 

(11.33) 

(T^^{t) ^ I dr [ dT'MG(i-r)iy(r,r')G^(t-T'), 
Jo Jo 

(11.34) 

<T^P(<)= [<T?^(<)]^ , (11.35) 

o-V{i)= [ dr [ dr' MG{t~T)u{T,T')G'^{t-T')M . 
Jo Jo 

(11.36) 

As explained above, we have G^{t,T) — G(i,T), so that 
we do not need to consider the transpose of the propaga- 
tor in Eqs. (II.33)-(II.36). 

Strictly speaking, since we already have the solutions, 
we do not need any master equation, but we provide it for 
completeness. The derivation of Ref. [1] can be staight- 
forwardly generalized to the case of N system oscillators 
and the resulting master equation can be written as 

(11.37) 

given the coefficient matrices 

n{t) = -i{t)^-\t), (11.38) 

B{t) = i {n{t) arit) + cTrit) H^{t) + &T{t)} , 

(11.39) 

where H are the homogeneous coefficients which contain 
system renormalization and dissipation, and D are the 



diffusive and anti-diffusive coefficients. (Note that the 
invalidity of previous derivations of the master equation 
for nonlocal dissipation pointed out in Ref. [1] applies in 
the same way to this case.) 



III. REFINEMENT OF THE SOLUTIONS 

A. Local Damping 

For local dissipation it is very easy to solve for the 
phase-space propagator, as it is formally determined by 
the matrix exponential 

*(i) =e-*«, 

Ml 

C 27oM-i 



n 



(III.l) 

(III.2) 



Directly solving for the position propagator involves an 
equivalent but quadratic matrix equation. Even for large 
systems, diagonalization of ft perturbatively in 7g is 
straightforward. 



B. Rational Damping 

1. Pseudo-normal Modes 

If one has a rational damping kernel in the Laplace 
domain (or even meromorphic) , then the propagator will 
also be rational (or meromorphic) in the Laplace domain. 
Therefore, it can be decomposed into simple fractions 
which correspond to pseudo-normal modes: 



G(s) = 
G{t) = 



(III.3) 
(III.4) 



assuming no repeated roots, and where the Ufe form an 
overcomplete basis given the initial condition G(0) = 
M~i, symmetry and the uniqueness of solutions. These 
are not true modes and the Ufe are not orthonormal be- 
cause for nonlocal damping the pseudo-modes outnumber 
the system degrees of freedom. However, they will be as- 
sociated with non-physical modes in a higher-dimensional 
linear system in the following subsection. 

From the frequency representation of Eq. (11.18) and 
the Langevin equation, the pseudo-mode decomposition 
must satisfy 



2sA(s) -I- n- 



s- fk 



1 , (III.5) 



for all s. Taking the limit s ^ f^. reveals a necessary 
condition for convergence, the characteristic equation: 



U = 0, 



(III.6) 
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which is a nonlinear eigen-value equation. Canonical- 
Uke perturbation theory apphed to this equation will be 
in agreement with results from the second-order master 
equation [m], but this method is much more efficient in 
calculating these particular frequencies to higher order 
as this method does not require any integration. 



2. Rational Damping: Extended Phase Space 

For the nonlocal damping kernel which is rational in 
Laplace space, we can transform the nonlocal problem 
into an effectively time-local and time-homogeneous one 
by considering a higher dimensional phase space. This 
was possible for one system oscillator, but is much more 
useful with N oscillators. The specific problem that we 
will work out will be the regulated Ohmic coupling 



7o 
1 + i 



(III.7) 



where the Ohmic-limiting damping matrix 7o ^^'^ f^e 
cutoff matrix A commute and are both positive definite. 
This damping kernel is a rational function of order [0/1] 
in s and corresponds to 



7(-) = 



(Ill.i 



in the Fourier domain; see Eq. (111.22) below. For this 
to be a valid coupling, one only needs to ensure that 
jiuj) > for all w, which is indeed the case. 

Together with the initial conditions (Xo,Po) 
Eq. (11.11) yields the following initial-value constraints: 



X(0) = Xo , 
MX(0) - Po , 

MX(0) = - [C + 2A7o] Xo , 



(III.9) 
(111.10) 
(111.11) 



Q,a 



where the initial acceleration will always appear to have 
a bare or non-renormalized frequency. In general, if the 
damping kernel is of order [n/d] then we will need to 
determine d additional initial conditions beyond the first 
two. 

Next, we factor out the damping kernel's denomina- 
tor and represent the homogeneous Langevin equation in 
Laplace space as a polynomial 

[(1-l-sA-i) (s^M + C) +2.s7o] X(s) 

= (1 + sA-i) (sMXo + Po) , (111.12) 

which is equivalent to the third-order differential equa- 
tion 

:0, (IIL13) 



with constant coefficients. Given that we have made no 
approximations, solutions to this equation with proper 
initial conditions are immune to any issues of runaway 
solutions or causality violation often associated with 
higher-order differential equations. 

Because we have a linear system of ordinary differential 
equations, we can solve this system with a matrix expo- 
nential analogous to the case of local damping. Here it 
is natural to consider the extension of phase space into 
"forces" such that 



d 
dt 



X 
P 
F 








M-1 




X" 


1 




P 


2A7o-C)M-i -A 




F 



(III.14) 

and let us denote this effective time-translation generator 
as so that we can write our solutions as 



(III.15) 



where the initial condition Fo is determined by 
Eq. (111.11). Given our initial value constraints, we can 
then express our solutions as 



-xit)- 




"Xo" 


P{t) 




Po 


m 




Fo 



X(t) 




Pit) 




m 





1 



-2A70-C 



Xo 

Po 




(III.16) 

The initial- value constraint matrix, which will be denoted 
by T, maps the two initial conditions (Xo, Pq, 0, • • • ) 
of the nonlocal equation into the larger number of ini- 
tial conditions required by the effective local equation 
(Xo, Pq, Fo, • • • ).^ Given this representation we can 
now identify the 2x2 leading principal minor of e* T 
with 4*(t), since this matrix would map (Xo,Po) onto 
(X(t),P(i)). In the particular example that we are con- 
sidering here Fq has no dependence on Pq and the second 

column of e*^ directly yields (G(t), MG(t), MG(i)) , 

so that one can read the result for the propagator and 
its first two derivatives and obtain $(i) making use of 
Eq. (11.28). Hence, in the end one has been able to find 
the matrix propagator for the integro-differential equa- 
tion by solving a higher-order linear differential equation. 



C. Late-Time Thermal Covariance 

Using the results in Appendix E of Ref. [1], we can 
reduce the late-time covariance to a single frequency in- 



^ In general, the dimension of the "extended phase space" and the 
ranlc of T will correspond to the order of the equivalent ordinary 
differential equation. 
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tegral in the Laplace domain: 



+00 



1 1 / ■ ' . ~ 

0-7^(5) w / duj ^{s+iuj)'N{uj) ^ {s — iuj) 

s 2tt , 



which is exact in the late-time limit 



(Tt{oo) — lim s a-T{s) 

s— >o 



(III.17) 



(III.18) 



and for all times if the dissipation is local. 

Moreover, there is a simplification to this expression 
which was not exceptionally useful for a single system 
mode but is extremely useful for multiple system modes. 
The expression can be further simplified from a quadratic 
form into a linear form. First we express the noise kernel 
in terms of the damping kernel by means of the FDR 
(11.12). Next we relate the Fourier transform of the 
damping kernel to its Laplace transform by first noting 



7(0;) 



+00 



dt e 



dt e 



--fit) 



dt e 



-\-iujt 



(III.19) 



7(-i). (in.20) 



Then we use the fact that the damping kernel is double- 
Hermitian (see Ref. [9]) to obtain 



7(^)= / die— *7(t) 
Jo 



dte 



Yit)- (111.21) 



Assuming no singularities in 7(s) for Re[s] > 0, which 
should be the case as the damping kernel is a somewhat 
localized distribution in time, we may analytically con- 
tinue the above relation to obtain 



7(w) = lim 7(s + iw) +7 [s — iu]) 
s— >o 



(III.22) 



Then we expand Re[7(za;)] pairing 'y[+iuj) with $(-|-iw) 

and 7 (— «w) with $ {—v-^)- One can then apply the 
Langevin equation evaluated at imaginary frequencies to 
eliminate any explicit dependence on the damping kernel, 
so that only the implicit dependence of the propagator 
is left. After some simplification and application of sym- 
metries we are left with the expression 



27r ./ ^ LU 



G{-~iuj) 

w2MG(-ia;)M 

(III.23) 



which avoids all phase-space matrix multiplication. 



IV. SYSTEM MODES IN A COMMON 
ENVIRONMENT 

Here we consider identical Unruh-DeWitt detectors Qi 
in a shared environment (a massless scalar field in a vac- 
uum or thermal state), but at different locations as 



considered in Ref. ['']. The coupling is local and linear 
between the system positions (more typically momenta) 
and field operator. In this model the damping kernel is 
given by 



7ij(iu) = 27osinc(njw) . 
%{s) = 70 ■ 



ras 



(IV.l) 
(IV.2) 



in the Fourier and Laplace domains respectively, and 
where r^j — Vi~Yj. The Markovian limit is reached here 
when the environment is at high temperature and the de- 
tectors are taken to be very far apart, (in which case one 
recovers 70 as the damping rate of the individual system 
oscillators) or very close together (in which case one has 
effectively local cross-damping) . At finite temperature it 
is well known that QBM requires a finite cutoff regula- 
tor for Ohmic coupling to the environment, and this is 
also the case here. Note that distance variation natu- 
rally takes the form of an Ohmic regulator in the limit of 
small separations. Following Ref [ ], which considered 
the same environment, we will not evaluate the correla- 
tions more precisely than to some smallest scale ro in the 
sense that we take 



lim 



A- 



(IV.3) 



to provide a natural cutoff of A. The scale tq can be 
interpreted as the effective size of the detectors. A more 
detailed analysis would require including the form factors 
for the detectors, but the conclusions should be qualita- 
tively the same. 

For a pair of identical detectors Qi and Q2, the dy- 
namics is naturally described in terms of the sum and 
difference Q± = Qi ± Q2 (which correspond to the nor- 
mal modes of the free theory) and a similar combination 
for the components of the damping kernels. 



7±(s) = 7ii(s) ± 712(5) 



(IV.4) 



as well as the noise kernel (by the FDR). Hence, when 
the two detectors are very close, their difference Q- ex- 
periences vanishing dissipation. It is associated with 
a decoherence-free subspace, involving dark states, but 
only in the limit of minimal separation so that their en- 
vironments appear identical and are perfectly correlated. 
This is essentially the regime considered in Ref. [ ] for the 
entanglement dynamics of two oscillators in a common 
environment. But more generally, for small separations 
the sum (5+ relaxes quickly and then on a much longer 
timescale the difference Q_ will also thermalize. The re- 
sultant asymptotic state, which is Gaussian in the sum 
and difference, can also be entangled at higher order in 
the coupling. 

More generally Ref. [ ,] considered finite separation dis- 
tances and determined the oscillators to be asymptoti- 
cally entangled when near the cutoff scale in separation. 
Our more general formalism allows for non-resonant de- 
tectors which cannot be transformed into a pair of indi- 
vidual quantum Brownian oscillators. Hence, it enables 
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us to determine how resonant the detectors must be for 
this asymptotic entanglement to ensue. 

A. Regulation and Integration 

The non-perturbative late-time covariance (III. 23) 
evaluated with our exponential cutoff-regulator (IV. 2), 
is exceptionally oscillatory and does not lend itself to 
well-behaved numerics, especially in the near regime. In 
evaluating the late-time covariance, the exact regulator 
in the Laplace domain 

X{z) = , (IV.5) 

z 

must be evaluated at imaginary arguments, and it be- 
comes both oscillatory and asymptotically 0{l/iz) as 
\z\ — oo. Therefore, the appropriate small-r approxima- 
tions of this regulator are the X[n/n+i] {z) Pade approxi- 
mants. These are the best local (small-z expansion) and 
rational (ratio of polynomials) approximations which are 
also 0{\/iz). We plot the first three such approximants 
in Fig. 1 . One can see that they are good approximations 



Re[Ar('z)l 




FIG. 1. ■ ■ ■ The first three X[n/n+\\{z) Pade approximants of 
the • exact regulator. 

up to r ^ n7r/2i7, similar to the Taylor series convergence 
for sinusoidal functions. In this work we have primarily 
considered the first-order Pade approximant 



X[o/i]{z) = 




(IV. 6) 



which is also an Ohmic regulator for small separations. 
This local approximation is most accurate in the high- 
cutoff and small-r regime, which is precisely where we 
need to investigate more carefully, but it also has the 
correct asymptotics. Resultant calculations can be ob- 
tained exactly and are very well behaved. In addition to 
numeric solutions and analytic small-r solutions, we also 
employ analytic small-70 solutions using weak-coupling 
perturbation of Eq. (III. 6). 



B. Entanglement Dynamics 

In calculations of entanglement we consider both the 
Peres-Horodecki criterion E [11] and the logarithmic neg- 
ativity Ej^ [12] as a consistency check. In all graphs, 
positive values denote entanglement and negative ones 
separability. We consider the zero-temperature regime 
only as it emphasizes entanglement. 

Our generic analysis of entanglement dynamics in- 
volves two factors: (1) the relevant timescales for deco- 
herence and (2) unmaximized entanglement monotones 
of the asymptotic state [1^]. In this linear model de- 
coherence arises mainly due to the growth of the ther- 
mal covariance ctt, which smears away oscillations in the 
Wigner function corresponding to quantum interference. 
A more thorough explanation of the general evolution is 
given in Sec. 4 of Ref. [±]. 

En 
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FIG. 2. Qualitative plot of (unmaximized) entanglement dy- 
namics including entanglement sudden death, revival, and 
asymptotic separability. 

Traditionally defined entanglement monotones are not 
sufficient to give a cursory analysis of entanglement sud- 
den death as they do not distinguish between separable 
states. Given some operationally defined entanglement 
monotone £ = max[0,£], so that £ is zero for all sepa- 
rable states, then the unmaximized function £ is more 
useful when its evolution is continuous. As we quali- 
tatively plot in Fig. 2, entanglement sudden death oc- 
curs because the unmaximized entanglement monotone 
asymptotes towards a negative value, whereas the tradi- 
tionally defined entanglement monotone does not asymp- 
tote towards zero. Therefore, given that decoherence in 
linear models is fairly well understood, we will focus pri- 
marily upon analysis of the asymptotic (unmaximized) 
entanglement monotone. 



1. Close Detectors 

In agreement with Ref. [''] we find the resonant oscilla- 
tors to be asymptotically entangled when their separation 
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FIG. 3. The asymptotic (unmaximized) entanglement mono- 
tones of two resonant oscillators as a function of separation 
where 70 = Q/W and A = WOQ. 



is comparable to the inverse of the cutoff scale. Fig. 3 de- 
notes the asymptotic entanglement measures of the reso- 
nant oscillators on the near cutoff scale where asymptotic 
entanglement can emerge. It is curious that, despite the 
same environmental correlations and despite the emer- 
gence of a dark state under the same conditions, for a 
pair of two-level atoms the possibility of asymptotic en- 
tanglement does not emerge (at second-order) for perfect 
cross-correlations [10]. Moreover the behavior is rather 
opposite to that here, with the unmaximized entangle- 
ment monotones becoming more negative at proximity. 



2. Off-Resonant Detectors 

We know that asymptotic entanglement will ensue 
when the two oscillators are very close together and at 
resonance. Therefore the question arises, how close and 
how at resonance must the two oscillators be to have 
asymptotic entanglement. The first question has been 
answered: to have asymptotic entanglement, the oscilla- 
tors must be extremely close together, of the order of the 
inverse cutoff. The second question requires a more gen- 
eral multivariate treatment, which we can provide with 
this formalism. 

Here we take the two oscillators to be very close in 
position but with different frequencies fl± = Qq ± SQ. 
Expressed in this manner with average frequency Qq and 
difference 2(517, we find the asymptotic entanglement to 
be very much insensitive to SQ. Fig. 4 denotes the asymp- 
totic entanglement measures of the two oscillators for a 
moderate coupling strength. The stronger the environ- 
mental interaction is, the more asymptotic entanglement 
ensues and the less resonance sensitivity there is. But 
even for weak interactions, resonance does not appear 
to be a very stringent criterion for asymptotic entangle- 
ment, at least far less stringent than proximity. 
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FIG. 4. The asymptotic entanglement measures of two os- 
cillators as a function of detuning where 70 ~ O./10 and 
A = 100 



C. Sub- and Super-radiance 

We know that the dark and bright states emerge when 
the two oscillators are very close together and at res- 
onance. Therefore, the following question arises: how 
close and how close to resonance must the two oscillators 
be to have sub- and super-radiance? 



r 




10 20 30 40 

FIG. 5. Phenomenological decay rates as a function of sep- 
aration for two resonant oscillators using ■ weak-coupling 
perturbation and ■ ■ ■ the small-r Fade approximation where 

70 = n/10. 

In Fig. 5 we plot the phenomenological damping rates 
r as a function of separation between two resonant os- 
cillators. The sub-radiant decay rate starts at zero while 
the super-radiant decay rate starts doubled. Both then 
asymptote to the same value for large separations. For 
the exact regulator, the intermediate behavior of the 
decay rates is somewhat oscillatory and the separation 
timescales play a fairly different (and more causal) role 
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than a simple decay rate. But the important information 
to extract is that the existence of the dark and bright 
states is not particularly sensitive to a lack of proximity 
(where our approximation is valid). The relevant length 
scale appears to simply be the resonant wavelength. 

In Fig. 6 we plot the phenomenological damping and 
oscillation rates T and fZp as a function of detuning 
between two close oscillators' free parameters. To be 
more specific, these oscillators have effective mo des which 
evolve with frequencies F ± i fir where Clr = i/fip — F^. 
As with two-level atoms [10], the tuning needs to be bet- 
ter than the damping rate. Here we can also see that the 
frequency of the damped oscillators initially resists any 
display of detuning until this critical threshold is reached. 
Although the frequency fir is very similar for the bright 
and dark modes, they oscillate at different rates since 
their damping rates (and hence the frequency f2r) are 
very different. 

V. DISCUSSION 

In this work we have considered the general model of 
a system of quantum oscillators bilinearly coupled to an 
environment of quantum oscillators. We have shown how 
the phase-space solutions of Ref. [ ] for a single oscillator 
can be easily extended to this more general case. More- 
over, we have provided useful relations for the multipar- 
tite case, which greatly reduce the amount of necessary 
calculation. 



We have applied our formalism to the problem of local 
oscillators interacting with a common scalar-field envi- 
ronment. Our formalism is completely general and it has 
allowed us to explore this problem more thoroughly than 
before [ , ]. In particular we have found that detuning 
of the atoms does not have a profound effect upon the 
asymptotic entanglement. 

In Ref. [ ] we were able to provide exact, analytical 
solutions for various power-law spectral densities: sub- 
Ohmic, Ohmic, and supra-Ohmic. Multivariate envi- 
ronment correlations arising from a quantum field are 
considerably more complex, because the spatial sepa- 
ration between the local oscillators causes a functional 
dependence for the generalized damping kernel which 
is far more complicated than an approximately power- 
law distribution. In this work we have applied a three- 
way approach of (1) numerical integration of our simpli- 
fied expressions, (2) weak-coupling perturbation and (3) 
an appropriate Fade approximation of the influence ker- 
nel. Future work must push the boundaries of these ap- 
proximation schemes and consider alternatives, such as 
the perturbative cross-correlation expansion considered 
in Ref. [ ]. 
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FIG. 6. 



Homogeneous timescales as a function of detuning for two close oscillators where 70 



= n/100. 



